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1  EXECUTIVE  SUMMARY 


Several  computer  codes,  some  of  them  available  commercially  and  others  in  the  public 
domain,  currently  exist  which  can  predict  the  electronic  structure  of  molecules  of  interest. 
Similarly,  several  codes  exist  which  can  predict  the  kinetics  and/or  dynamics  of  a  system 
given  the  electronic  structure  of  the  molecules  involved.  Several  interfaces  between  some 
of  these  codes  currently  exist,  however,  none  are  seamless.  They  do  not  predict  the  likely 
error  based  on  the  level  of  quantum  theory,  nor  do  they  account  for  the  importance  of 
small  wells  in  the  reaction  entrance  or  exit  channels. 

SARA,  Inc.  has  teamed  with  the  University  of  Minnesota  (U  of  Minn.),  which  has 
developed  a  means  to  minimize  this  disconnect,  with  an  innovative  and  robust  design. 
Based  on  the  POLYRATE  modifications  developed  by  the  computational  chemistry 
group  at  the  U  of  Mirm.  led  by  Professor  Donald  G.  Truhlar,  the  redesigned  POLYRATE 
program  allows  users  to  accurately  predict  the  behavior  of  a  chemically  reacting  system. 

SARA,  Inc.  and  the  Department  of  Chemistry  and  Supercomputing  Institute  of  the 
University  of  Minnesota  -  Minneapolis  are  proud  to  present  the  results  of  the  Phase  1 
efforts  on  computational  predictions  of  kinetic  rate  constants.  The  purpose  of  these 
efforts  were  to  develop  seamless,  easy  to  use,  efficient  code  to  calculate  electronic  wave 
functions  and  potential  energy  surfaces  of  molecules  and  to  predict  kinetic  rate  constants 
for  reactions  a  priori.  From  these  modifications  to  POLYRATE,  we  have  demonstrated 
the  capability  to  calculate  kinetic  rate  constants  for  chemical  reactions  from  inputs 
obtained  by  electronic  structure  theory  with  seamless  integration  of  electronic  structure 
and  kinetics  calculations  as  described  in  this  final  technical  report  and  the  accompanying 
technical  papers. 

2  TECHNICAL  OBJECTIVES 

Work  on  the  project  started  in  August  2005.  Meetings  took  place  between  personnel  from 
SARA,  Inc.  (Drs.  Edward  Patton  and  Dmitry  Altshuler)  and  personnel  fi-om  the 
University  of  Minnesota  (U  of  Minn)  (Drs.  Don  Truhlar  and  Stephen  Klippenstein)  with 
the  kickoff  meeting  taking  place  on  September  8-9,  2005.  At  this  meeting,  also  attended 
by  several  U  of  Minn  post  doctoral  associates,  directions  for  the  project  and  roles  of 
participants  were  further  clarified  and  remaining  issues  pertaining  to  the  intellectual 
property  agreement  were  resolved.  Specific  goals  were  set: 

•  to  incorporate  new  kinds  of  dividing  surfaces  for  barrierless  association  reactions 
into  POLYRATE 

•  to  develop  improved  error  estimates  for  the  reliability  of  the  calculations 

•  to  improve  the  treatment  of  anharmonicity 

•  to  demonstrate  the  capability  for  convenient,  efficient,  and  reliable  direct 
dynamics 

•  to  make  generally  useful  improvements  in  POLYRATE  and  the  methods  used 
with  POLYRATE 

•  to  make  the  code  portable  under  Linux 


Contract  No.  FA9550-05-C-0065 


1 


It  was  determined  that  the  project  would  benefit  from  the  work  done  by  Dr.  Klippenstein 
who  had  been  working  with  Professor  Truhlar  on  a  related  project  on  POLYRATE. 
Because  of  this,  some  of  the  subtasks  initially  planned  were  deemed  not  necessary,  with 
the  previous  results  being  leveraged  into  this  effort.  Specifically,  the  original  intent  of 
these  efforts  had  been  to  enhance  POLYRATE  in  order  to  meet  two  of  the  three  technical 
Phase  I  objectives: 

•  to  include  error  estimates 

•  to  add  handling  of  certain  additional  types  of  chemical  reactions. 

When  the  proposal  was  written  originally,  it  appeared  that  additional  mathematical 
analysis  would  be  needed. 

But  based  on  the  previous  work  already  performed  by  Dr.  Klippenstein,  this  turned  out 
not  be  the  case.  In  the  original  plan  there  were  four  fundamental  subtasks  required  to  be 
accomplished  before  POLYRATE  would  be  ready  for  demonstration:  (1)  preliminary 
design,  (2)  handling  of  barrierless  reactions  and  reactions  with  potential  wells  must  be 
added,  (3)  error  estimates  were  to  be  added  as  a  result  of  mathematical  analysis,  and  (4) 
new  equations  were  to  be  implemented  in  POLYRATE.  It  was  determined  that  the 
subtasks  (2)  and  (3)  could,  in  fact,  be  incorporated  into  subtask  (4)  which  it  was  at  a 
savings  in  time  and  resources. 

The  preliminary  design  of  the  software  modifications  came  out  of  that  same  meeting  held 
on  September  8-9,  2005  at  the  University  of  Minnesota.  In  addition,  a  Master  Equation 
solver  was  determined  at  that  time  as  feasible  to  be  included.  These  enhancements  would 
enable  POLYRATE  to  handle  reactions,  other  than  the  simple  barrier.  It  was  also 
determined  that  the  U  of  Minn,  investigators  already  had  good  handle  on  errors  that  occur 
in  such  calculations.  The  main  source  of  such  errors  is  in  the  calculation  of  potential 
energy  surfaces.  Estimates  for  these  errors  would  be  incorporated  into  POLYRATE. 

2.1  Approach 

The  original  POLYRATE  software  uses  some  modem  versions  of  transition  state  theory 
(TST),  such  as  canonical  variational  theory  (CVT),  improved  canonical  variational  theory 
(ICVT),  and  microcanonical  variational  theory  (// VT).  These  theories  are  described  in 
journal  articles  and  in  a  chapter  of  the  upcoming  volume  in  the  Comprehensive  Chemical 
Kinetics  series,  but  are  not  widely  known  by  non-experts.  Enhancements  in  POLYRATE, 
implemented  within  the  Phase  I,  include  solving  the  Master  Equation  and  treatment  of 
barrierless  reactions. 

Prediction  of  reaction  rates  in  chemical  systems  when  the  concentrations  of  the  reactants 
vary  over  time  is  a  difficult  and  challenging  problem.  While  codes  for  this  prediction 
have  been  developed,  they  are  exceedingly  complex.  Yet,  they  require  expertise  at  a  level 
where  they  cannot  be  used  by  designers  and  engineers  working  on  chemical  or 
combustion  processes.  In  addition,  these  codes  do  not  provide  good  estimates  of  error,  so 
even  if  a  designer  or  engineer  could  use  them,  they  would  not  know  how  accurate  their 
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predictions  were.  SARA,  teamed  with  the  University  of  Minnesota  (U  of  Minn.)  has 
provided  an  innovative  and  robust  solution  to  this  problem  based  on  the  POLYRATE 
code  developed  by  the  computational  chemistry  group  at  the  U  of  Minn.  The  solution 
presented  here  is  an  enhancement  of  the  POLYRATE  program  as  needed  to  meet  the 
stated  requirements. 

The  intent  of  the  project  was  to  develop  a  computational  tool  usable  by  a  chemist  or  a 
chemical  engineer  who  is  not  an  expert  in  chemical  kinetics  or  quantum  chemistry.  It 
was,  however,  reasonable  to  assume  that  the  potential  user  of  the  product  had  a  solid  one- 
year  undergraduate  course  in  physical  chemistry.  The  user  should,  therefore,  be  familiar 
with  the  basic  concepts  and  terminology  of  quantum  chemistry,  statistical  mechanics  and 
chemical  kinetics.  It  would  still  be  helpful  for  the  user  to  have  had  a  more  mathematically 
sophisticated  physical  chemistry  course  and/or  a  modem  course  in  chemical  kinetics  (e.g. 
based  on  the  textbook  by  Steinfeld,  Francisco,  and  Hase'),  but  it  should  not  be  a 
requirement. 

2.2  Statement  of  the  Problem 

To  review,  in  order  to  accurately  predict  the  behavior  of  a  chemically  reacting  system,  it 
is  critical  to  know  how  the  reaction  rate  varies  with  the  concentration  (or  partial  pressure 
for  gas-phase  reactions)  of  the  reactants.  This,  in  turn,  requires  knowledge  of  a  parameter 
called  the  reaction  rate  coefficient,  which  itself  depends  on  temperature.  Mathematically, 
the  reaction  rate  coefficient,  k,  may  be  approximated  by  the  well  known  Arrhenius 
equation: 


k  =  Ae-^"'^ 
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where  A  is  referred  to  as  a  pre-exponential  factor,  E  is  the  activation  energy,  E  is  the  gas 
constant,  and  T  is  temperature.  Therefore,  in  order  to  have  an  approximation  to  the  rate 
constant  over  a  limited  temperature  range,  calculate  the  reaction  rate  coefficient,  it  is,  in 
fact,  necessary  to  calculate  two  constants:  A  and  E,  hereafter  referred  to  as  Arrhenius 
parameters.  In  order  to  calculate  a  rate  constant  over  a  broad  temperature  range,  one 
needs  three  or  more  parameters. 

However,  in  order  to  calculate  a  rate  constant  over  a  broad  temperature  range,  one  needs 
a  more  complex  model  with  at  least  three  parameters,  or  possibly  more.  The  team 
members  at  the  University  of  Minnesota  worked  to  develop  a  complex  input  model  to 
address  these  additional  parameters.  This  work  was  largely  based  on  the  existing 
POLYRATE  program  originally  developed  by  U.  of  Minn.  Phase  I  work  has  focused  on 
extending  the  applicability  of  POLYRATE  to  a  wider  type  of  chemical  reactions  and  on 
implementing  error  estimates  in  the  calculations.  Preliminary  designs  for  the 
POLYRATE  software  interface  determined  that  some  previously  developed  code  could 
be  incorporated  directly  in  POLYRATE. 


'  Steinfeld,  J.I.,  Francisco,  J.S.,  and  Hase,  W.L.  Chemical  Kinetics  and  Dynamics,  2"''  ed.,  Prentice  Hall, 
1999. 
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The  most  successful  method  for  calculating  rate  constants  goes  back  to  the  work  of  H. 
Eyring  and  his  co-workers^  and  is  known  as  Transition  State  Theory  (TST),  which  has 
become  a  workhorse  of  chemical  kinetics.  TST  provides  mathematical  equations  for 
calculating  rate  constants  and  hence  Arrhenius  parameters.  It  is  based  on  the  structure  of 
the  so-called  transition  state  between  the  reactants  and  the  products,  which,  in  theory,  can 
be  predicted  from  the  structure  of  the  reactants.  For  simple  molecules  and  for  the  simplest 
version  of  the  theory,  as  used  by  Eyring,  these  calculations  on  a  qualitative  basis  can  be 
carried  out  analytically  by  hand. 

However,  more  complex  cases  and  the  full  modem  theory  (to  be  called  generalized  TST 
or  GTST),  which  is  much  more  accurate,  require  numerical  calculations  with  a  computer. 
Kinetic  calculations  must,  essentially,  accomplish  two  tasks:  Predict  the  stmcture  of  the 
transition  state  or,  more  generally,  the  part  of  the  potential  energy  surface  needed  for  the 
d5mamics  calculations  and  (2)  calculate  the  rate  constant  itself,  or,  if  one  prefers  to  state 
the  results  that  way,  calculate  the  Arrhenius  parameters,  from  the  GTST  equations.  Both 
subjects  are  active  fields  of  research,  but  they  have  progressed  to  the  point  where  useful 
computer  codes  exist  to  calculate  electronic  stmcture  of  molecules,  and  useful  computer 
codes  which  calculate  rate  constants  or  Arrhenius  parameters  from  the  known  electronic 
stmcture  of  the  molecules  involved  in  the  reaction  also  exist.  There  are  interfaces 
between  these  two  types  of  codes,  but  none  are  seamless,  and  the  codes  are  not  as  user 
friendly  as  one  would  like. 

Another  drawback  in  the  existing  computational  codes  is  lack  of  error  estimates.  In  order 
to  select  the  best  computational  method,  it  is  necessary  to  know  the  error  it  may 
introduce.  Furthermore,  existing  computational  methods  work  better  for  simple  barrier 
reactions  than  for  reactions  with  chemically  significant  potential  wells.  To  summarize, 
the  following  requirements  should  be  met  in  order  to  make  such  computational  codes 
more  useful  for  chemical  engineers: 


•  It  should  be  possible  to  enter  the  formulas  for  the  reacting  molecules  and  obtain 
reaction  rates  without  having  to  know  the  electronic  stmcture  of  the  molecules 
involved.  If  several  reactions  can  occur  between  these  reactants,  it  should  be 
possible  to  predict  rates  for  all  of  them. 

•  Since  several  methods  are  available  for  kinetics  calculations,  it  is  necessary  to 
estimate  an  error  involved  in  using  each  of  them. 


This  work  focused,  therefore,  on  improvements  to  the  usefulness  of  one  of  the  available 
computational  kinetics  programs,  POLYRATE,  developed  by  Professor  Tmhlar’s  group 
at  the  University  of  Minnesota. 


^  Glasstone,  S.,  Laidler,  K.J.,  and  Eyring,  H.  Theory  of  Rate  Processes,  New  York:  McGraw-Hill,  1941. 
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2.3 


POLYRATE  Program 


The  POLYRATE^  program  has  proven  to  be  very  successful  in  calculation  of  the 
Arrhenius  parameters  for  a  large  class  of  chemical  reactions.  POLYRATE  excels  for  the 
so-called  simple  barrier  reactions'*,  for  which  along  with  conventional  TST  it  uses  several 
varieties  of  the  Generalized  Transition  State  Theory  (GTST)^,  also  developed  by  Truhlar 
and  his  co-workers:  Canonical  Variational  Theory  (CVT),  Improved  Canonical 
Variational  Theory  (ICVT),  and  Microcanonical  Variational  Theory  (//VT).  These  are  all 
forms  of  Variational  Transition  State  Theory^  (VTST).  There  is  a  considerable  amount  of 
literature  on  these  methods.  Therefore,  they  are  not  described  in  depth  here. 

All  varieties  of  the  TST  postulate  the  existence  of  the  transition  state  (sometimes  also 
called  the  activated  complex),  which  acts  as  a  potential  barrier  between  the  reactants  and 
products.  Arrhenius  parameters  in  POLYRATE  are  calculated  by  finding  the  minimum 
energy  path  (MEP)  across  this  barrier.  As  the  name  implies,  this  postulate  is  certainly  true 
for  simple  barrier  reactions.  By  their  very  nature,  the  computational  methods  in 
POLYRATE  rely  on  the  electronic  structure  of  the  molecules  involved.  In  fact, 
POLYRATE  can  accept  as  input  an  electronic  structure  file.  Therefore,  at  least  for  simple 
barrier  reactions,  this  program  already  fulfils  the  first  of  the  stated  Phase  I  objectives  -  to 
demonstrate  the  capability  to  calculate  kinetic  rate  constants  for  chemical  reactions  from 
inputs  obtained  by  electronic  structure  theory.  However,  input  files  in  POLYRATE  must 
be  built  manually,  which  requires  a  considerable  amount  of  labor. 

The  Phase  1  work  focused,  on  one  hand,  on  improving  POLYRATE  for  reactions  other 
than  those  of  the  simple  barrier  type,  and,  on  the  other  hand,  on  estimating  errors  that 
may  be  present  in  computational  algorithms.  Phase  2  work  will  focus  on  eliminating  the 
manual  building  of  input  files  simplifying  the  interface  between  electronic  structure 
theory  and  dynamics  by  better  integrating  POLYRATE  with  electronic  structure  tools 
developed  by  Professor  Truhlar’s  group. 

2.4  Enhancements  to  POLYRATE 

Not  all  chemical  reactions  are  of  simple-barrier  type.  A  common  occurrence  is  that  one  or 
more  potential  wells  lie  along  the  reaction  path.  Should  the  system  fall  into  such  a  well,  it 
may  become  stabilized  instead  of  decomposing  into  the  reaction  products  or  back  into  the 
reactants.  Furthermore,  it  is  possible  for  a  system  in  such  a  well  to  isomerize  into  a 
system  in  another  well  (perhaps  more  than  once)  before  decomposing  into  reactants  or 


^  Description  of  POLYRATE,  including  the  complete  manual,  can  be  found  at 
http://comp.chem.umn.edu/nolvrate/ 

^  Femandez-Ramos,  A.,  Miller,  J.A.,  Klippenstein,  S.J.,  and  Truhlar,  D.G.  “Bimolecular  Reactions”,  in 
Comprehensive  Chemical  Kinetics,  R.W.  Carr  and  M.R.  Zachariah,  eds.,  to  be  published.  Simple  barrier 
reactions  and  computational  algorithms  for  them  are  described  in  Section  2.4.  Other  types  of  reactions  are 
treated  in  Section  2.5. 

^  Truhlar,  D.G.,  Isaacson,  A.D.,  and  Garrett,  B.C.  “Generalized  Transition  State  Theory”  in  The  Theory  of 
Chemical  Reaction  Dynamics,  M.  Baer,  ed.,  Boca  Raton:  CRC  Press.,  1985,  Vol.4,  pp.65-137. 

®  Truhlar,  D.G.  and  Garrett,  B.C.  “Variational  Transition  State  Theory”  Rev.  Phys.  Chem.,  1984, 
Vol.35,pp.l59-189. 
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products.  It  can  also  decompose  into  several  sets  or  reaction  products.  There  is  also  a 
possibility  that  a  potential  energy  barrier,  postulated  in  TST,  is,  in  fact,  absent.  Barrierless 
reactions  with  multiple  wells  are  a  key  class  of  reactions  in  practical  applications. 

2.4.1  Barrierless  Reactions 

VTST  is  applicable  even  to  barrierless  reactions'^.  However,  the  rigid-rotor  harmonic- 
oscillator  treatment  of  molecules  turns  out  to  be  very  inaccurate  for  transition  states  of 
such  reactions.  Therefore,  it  is  necessary  to  introduce  the  treatment  of  anharmonicity  into 
the  calculations.  This  requires  major  changes  including  new  coordinate  systems  to  define 
the  dividing  surfaces. 

At  the  same  time,  the  absence  of  a  barrier  allows  for  some  simplifications.  One  of  them  is 
that  quantum  tunneling  effects  are  not  very  significant.  Other  applicable  simplifications 
are  described  in  depth  in  the  above-cited  reference''.  Work  included  implementation  of 
these  calculations  in  POLYRATE  as  well  an  estimation  of  possible  errors. 

2.4.2  The  Master  Equation 

Chemical  reactions  involving  multiple  wells  or  multiple  sets  of  reaction  products  are 
studied  using  the  statistical  model  governed  by  the  master  equation.  In  its  most  primitive 
form  the  master  equation  is  written  as 

where  ni(t)  is  the  population  of  the  state  i  at  the  time  t,  and  pij  is  the  probability  of 
transition  from  state  j  to  state  i.  Essentially,  the  master  equation  describes  a  stochastic 
Markov  process  as  long  as  the  coefficients  py  do  not  depend  explicitly  on  time  or  the  past 
history  of  the  populations.  This  assumption  is  satisfied  for  all  problems  of  interest. 
Furthermore,  it  is  often  adequate  to  assume  that  these  coefficients  do  not  depend  on  ni(t), 
which  makes  the  equation  linear. 

Some  alternative  approaches  to  solving  the  Master  Equation  were  pursued  as  part  of  this 
project.  The  importance  of  Markov  processes  in  other  areas  of  engineering  and  operations 
research  led  to  the  development  of  several  algorithmic  approaches  for  their  analysis^. 
While  most  of  these  approaches  dealt  mainly  with  queuing  theory  and  decision  processes, 
the  underlying  mathematical  methods  may  be  applicable  to  any  continuous-time  Markov 
processes. 

In  addition,  some  alternative  approaches  to  solving  the  master  equation  were  be  pursued 
as  part  of  this  project.  The  importance  of  Markov  processes  in  other  areas  of  engineering 
and  operations  research  led  to  the  development  of  several  algorithmic  approaches  for 


^  See,  for  example,  Ross,  S.M.  Applied  Probability  Models  with  Optimization  Applications.  New  York: 
Dover,  1 992;  Neuts,  M.F.  Matrix-Geometric  Solutions  in  Stochastic  Models:  An  Algorithmic  Approach. 
New  York:  Dover,  1994;  Tijms,  H.C.  Stochastic  Models:  An  Algorithmic  Approach.  Chichester:  John 
Wiley  &  Sons,  1994. 
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their  analysis*.  While  most  of  these  approaches  deal  mainly  with  queuing  theory  and 
decision  processes,  the  underlying  mathematical  methods  may  be  applicable  to  any 
continuous-time  Markov  processes. 

2.4.3  Graphical  User  Interface  (GUI) 

Not  all  chemical  reactions  are  of  simple-barrier  type.  A  common  occurrence  is  that  one  or 
more  potential  wells  lie  along  the  reaction  path.  Should  the  system  fall  into  such  a  well,  it 
may  become  stabilized  instead  of  decomposing  into  the  reaction  products  or  back  into  the 
reactants.  Furthermore,  it  is  possible  for  a  system  in  such  a  well  to  isomerize  into  a 
system  in  another  well  (perhaps  more  than  once)  before  decomposing  into  reactants  or 
products.  It  can  also  decompose  into  several  sets  or  reaction  products.  As  discussed 
above,  there  is  also  a  possibility  that  a  potential  energy  barrier,  postulated  in  TST,  is,  in 
fact,  absent.  Barrierless  reactions  with  multiple  wells  are  a  key  class  of  reactions  in 
practical  applications. 

To  access  these  newer  levels  of  complexity  being  built  into  POLYRATE,  a  sequence  of 
GUI  screens  were  designed.  Each  screen  was  designed  to  eontain  necessary  menu  options 
with  the  available  new  choices.  Chemical  species  (both  reactants  and  products)  could  be 
entered  by  drawing  their  structures  with  a  mouse.  For  instance,  if  the  user  decides  to 
calculate  potential  energy  from  scratch,  the  calculation  could  be  accomplished  by  a 
seamless  interface  with  one  of  the  electronic  structure  packages  developed  by  the 
University  of  Minnesota,  e.g.,  GAMES SPLUSRATE,  which,  in  turn,  calls 
GAMESSPLUS,  another  package  developed  by  the  University  of  Minnesota. 

It  is  of  importance  to  note  that  no  specialized,  expert-level  knowledge  of  PES  is  needed 
to  provide  input  for  this  section.  Progress  was  also  made  on  making  POLYRATE 
portable,  to  run  under  LINUX,  specifically,  under  g77  compiler,  the  standard  compiler 
for  LINUX  platforms.  Based  on  analysis  and  assessment  of  electronic  structure  tools  also 
developed  by  the  same  group  at  the  University  of  Minnesota,  the  final  deliverable  of  the 
Phase  I  will  be  a  system-level  specification  of  the  new  integrated  tool.  Phase  II  work  will 
proeeed  to  actually  design,  test  and  implement  the  tool  from  this  specification. 

At  the  same  time,  it  was  decided  that  it  would  be  unreasonable  not  to  give  the  user  an 
option  to  select  the  computational  approach.  In  fact,  the  user  may  want  to  try  more  than 
one  such  approach  in  order  to  compare  the  results.  Furthermore,  while  POLYRATE  gives 
the  user  options  for  various  treatments  of  anharmonicity  and  tunneling,  it  was  decided 
that  it  was  unreasonable  to  take  these  options  away  for  the  sake  of  simplifying  the 
usability  of  the  produet.  Therefore,  it  seemed  that  the  best  option  was  to  allow  the  user  to 
select  a  computational  approach  from  the  available  choices,  but  to  include,  with  the 
software,  a  comprehensive  description  of  these  methods.  Such  a  description  should 
probably  take  the  form  of  an  advanced  textbook,  rather  than  just  a  software  reference 


*  See,  for  example,  Ross,  S.M.  Applied  Probability  Models  with  Optimization  Applications.  New  York: 
Dover,  1992;  Neuts,  M.F.  Matrix-Geometric  Solutions  in  Stochastic  Models:  An  Algorithmic  Approach. 
New  York:  Dover,  1994;  Tijms,  H.C.  Stochastic  Models:  An  Algorithmic  Approach.  Chichester:  John 
Wiley  &  Sons,  1994. 
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manual  and  would  be  aimed,  as  stated  earlier,  at  a  reader  with  a  solid  undergraduate 
coursework  in  physical  chemistry. 

Presently,  in  order  to  use  POLYRATE,  the  user  must  manually  build  an  input  file.  In  the 
final  Phase  2  form,  this  will  be  replaced  by  a  sequence  of  screens.  Each  screen  will 
contain  menu  options  with  available  choices.  Chemical  species  (both  reactants  and 
products)  will  be  entered  by  drawing  their  structures  with  a  mouse. 

3  RESULTS 

The  following  subsections  describe  how  the  POLYRATE  input  files  will  be  built,  based 
on  this  Phase  1  software  design.  The  modifications  to  the  algorithms  that  were  developed 
for  this  effort  and  the  references  to  the  papers  under  which  these  results  are  given  in 
detail.  Electronic  copies  of  all  papers  being  generated  and  citing  these  results  will  be 
included  on  CD-ROM  and  furnished  to  the  Government.  To  seamlessly  access  these 
modifications  as  part  of  the  entire  POLYRATE  program,  the  subsequent  sections  of  this 
report  detail  the  modified  GUI  designed  to  access  both  existing  and  new  modes  of  the 
software.  These  sections  are  written  from  the  internal  point  of  view,  i.e.,  they  are 
organized  around  the  sections  and  keywords  of  the  existing  POLYRATE  main  input  file. 

3.1  B  imolecular  reactions 

Initial  modifications  to  the  algorithms  dealt  with  the  theoretical  and  computational 
modeling  of  bimolecular  reactions,  especially  with  generally  applicable  methods  for 
kinetics  (i.e.,  overall  rates  as  opposed  to  detailed  dynamics).  These  included  a  basic 
theoretical  framework  that  could  be  used  for  gas-phase  thermal  reactions,  gas-phase 
microcanonical  and  state-selected  reactions,  and  condensed-phase  chemical  reactions. 

The  treatment  of  gas-phase  thermal  reactions  included  separate  discussions  of  simple 
direct  reactions  over  a  barrier,  which  usually  have  tight  transition  states  and  reactions 
proceeding  over  a  chemical  potential  well,  which  now  can  have  a  number  of  additional 
complications,  such  as  barrierless  addition  potentials  (which  generally  have  loose, 
flexible  transition  states),  competitive  reaction  pathways,  isomerizations  between 
multiple  wells,  and  pressure-dependent  energy  transfer  processes. 

With  a  heavy  emphasis  on  thermal  reactions  using  (generalized)  transition  state  theory 
(TST)  including  multidimensional  tunneling,  we  showed  that  this  provides  the  best 
available  method  to  calculate  thermal  rate  constants  for  all  but  the  very  simplest  systems. 
Work  on  state-selective  reactions  and  product  state  distributions  included  an  introduction 
to  the  theory  of  electronically  nonadiabatic  reactions  and  coupled  potential  energy 
surfaces,  as  required  for  modeling  photochemical  and  chemiluminescent  reactions.  For 
bimolecular  reactions  in  liquid  solution,  we  considered  diffusion  control  and  equilibrium 
and  nonequilibrium  solvation.  From  these  results,  there  has  been  great  progress  in  our 
ability  to  model  the  kinetics  of  bimolecular  reactions.  This  derives  from: 
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•  improved  methods  for  generating  and  using  reactive  potential  energy  surfaces, 
especially  implicit  potential  energy  surfaces  generated  by  direct  dynamics 

•  improved  dynamical  algorithms,  including  practical  methods  for  finding 
variational  transition  states,  well-validated  multidimensional  methods  for 
including  tunneling,  and  master  equation  methods  for  treating  nonequilibrium 
distributions,  especially  in  multi-well,  multi-arrangement  reactions 

•  efficient  methods  for  interfacing  both  of  the  above. 

We  anticipate  continued  improvements  in  the  theoretical  development  for  extending 
POLYRATE  to  further  barrierless  reactions  and  reactions  with  potential  wells^. 

3.2  Transition  states  with  wide-amplitude  motion 

Practical  approximation  schemes  for  calculating  partition  functions  of  torsional  modes 
were  tested  against  accurate  quantum  mechanical  results  for  H2O2  and  six  isotopically 
substituted  hydrogen  peroxides.  The  schemes  were  classified  on  the  basis  of  the  type  and 
amount  of  information  that  is  required.  First,  approximate  one-dimensional  hindered- 
rotator  partition  functions  were  benchmarked  against  exact  one-dimensional  torsion 
results  obtained  by  eigenvalue  summation.  The  approximate  one-dimensional  methods 
tested  in  this  stage  included  schemes  that  only  require  the  equilibrium  geometries  and 
frequencies,  schemes  that  also  required  the  barrier  heights  of  internal  rotation,  and 
schemes  that  required  the  whole  one-dimensional  torsional  potential.  Then,  three  classes 
of  approximate  full-dimensional  vibrational-rotational  partition  functions  were  calculated 
and  were  compared  with  the  accurate  full-dimensional  path  integral  partition  functions. 
These  three  classes  were: 

•  separable  approximations  combining  harmonic  oscillator-rigid  rotator  models 
with  the  one-dimensional  torsion  schemes 

•  almost-separable  approximations  in  which  the  nonseparable  zero-point  energy 
was  used  to  correct  the  separable  approximations 

•  improved  nonseparable  Pitzer-Gwinn-type  methods  in  which  approaches  of  type 
1  were  used  as  reference  methods  in  the  Pitzer-Gwinn  approach. 

The  effectiveness  of  these  methods  for  the  calculation  of  isotope  effects  was  studied. 
Based  on  the  results  of  these  studies,  the  best  schemes  of  each  type  were  recommended  in 
Reference  9  for  further  use  on  systems  where  a  corresponding  amount  of  information  is 
available. 

Neither  the  heights  nor  spacings  of  the  torsion  barriers  are  equal  for  the  torsional  motion 
in  H2O2  so  the  system  is  an  excellent  case  to  use  in  testing  methods  intended  for  general 
problems.  H2O2  is  also  the  only  system  with  a  torsion  of  any  kind  for  which  accurate, 
full-dimensional  quantal  partition  functions  have  been  published  for  use  as  benchmarks. 


®  "Modeling  the  Kinetics  of  Bimolecular  Reactions,"  A.  Femandez-Ramos,  J.  A.  Miller,  S.  J.  Klippenstein, 
and  D.  G.  Truhlar,  Chemical  Reviews  106, 4518-4584  (2006). 
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Twenty-two  methods  for  calculating  partition  functions  of  one-dimensional  torsional 
potentials  were  examined.  If  the  potential  can  be  easily  represented  in  a  low-order  Fourier 
series,  the  Hamiltonian  takes  the  form  of  a  banded  matrix  with  analytically  determinable 
matrix  elements;  for  many  applications  such  systems  can  often  be  affordably 
diagonalized  and  partition  functions  obtained  via  direct  eigenvalue  summation. 

When  this  approach  is  not  feasible,  but  one  can  still  afford  to  perform  a  numerical 
quadrature  over  the  torsion  coordinate,  the  Torsional  Displaced-Points  Path  Integral 
method-  with  Harmonic  Sampling  (TDPPI-HS)  is  the  method  of  choice.  If  the  potential  is 
only  known  at  the  locations  of  the  barriers  and  minima,  then  the  segmented  reference 
Pitzer-Gwinn  scheme  or  the  Segmented  Reference  potential  version  of  the  TDPPI-HS 
(SR-TDPPI-HS),  (which  requires  accurate  integration  involving  only  the  segmented 
reference  potential),  are  the  most  appropriate  methods  from  among  those  considered  here. 
For  systems  where  only  the  geometries  and  frequencies  are  known,  the  CT-Cco  method  is 
best  option.  We  also  point  out  that  the  usual  manner  of  applying  the  harmonic 
approximation  would  be  in  error  by  approximately  a  factor  of  2,  but  once  one  recognizes 
the  issue,  this  error  is  easily  corrected  without  requiring  additional  data. 

Zero-point  energy  (ZPE)  anharmonicity  has  a  large  effect  on  the  accuracy  of  approximate 
partition  function  estimates.  If  the  accurate  ZPE  is  taken  into  account,  separable 
approximation  partition  functions  using  the  most  accurate  torsion  treatment  and  harmonic 
treatments  for  the  remaining  degrees  of  freedom  agree  with  accurate  QM  partition 
functions  to  within  a  mean  accuracy  of  9%.  If  no  ZPE  anharmonicity  correction  is  used, 
the  mean  accuracy  of  the  separable  treatment  drops  to  23%.  We  also  presented  a  simple 
modification  of  the  multidimensional  Pitzer-Gwinn  approximation,  denoted  the 
improved-reference  Pitzer-Gwinn  (IRPG)  approximation,  that  includes  an  accurate 
treatment  of  a  ID  torsional  model  instead  of  a  harmonic  treatment  in  the  quantum 
correction  to  the  full-dimensional  classical  phase-space  integral.  When  an  accurate  ZPE 
correction  is  also  employed,  the  IRPG  method  provides  substantial  improvements 
compared  with  the  standard  Pitzer-Gwinn  approximation  with  the  worst  errors  for  H2O2 
being  only  1%'”. 

3.3  Direct  dynamic  calculations 

Further  modifications  included  an  efficient  method  for  interfacing  POLYRATE  with 
direct  dynamics  calculations,  which  is  an  important  element  in  making 
the  package  easier  to  use  by  nonspecialists.  Multiconfiguration  molecular  mechanics 
(MCMM)  is  a  general  algorithm  for  constructing  potential  energy  surfaces  for  reactive 
systems%' ' . 


"Statistical  Thermodynamics  of  Bond  Torsional  Modes.  Tests  of  Separable,  and  Almost-Separable,  and 
Pitzer-Gwinn  Approximations,"  B.  A.  Ellingson,,  V.  A.  Lynch,  B.  A.  Ellingson,  S.  L.  Mielke,  and  D.  G. 
Truhlar,  Journal  of  Chemical  Physics  125,  84305/1-17  (2006). 

"  Kim,  Y.;  Corchado,  J.  C.;  Villa,  J.;  Xing,  J.  and  Truhlar,  D.  G.  J.  Chem.  Phys.  2000,  112,  2718). 
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Our  results  show  how  the  performance  of  the  MCMM  method  can  be  improved  by 
adopting  improved  molecular  mechanics  parameters.  We  carried  out  calculations  of 
reaction  rate  constants  using  variational  transition  state  theory  with  optimized 
multidimensional  tunneling  on  the  MCMM  potential  energy  surfaces  for  three  hydrogen 
transfer  reactions,  and  we  compared  the  results  to  direct  dynamics.  We  found  that  the 
MCMM  method  with  as  little  as  one  electronic  structure  Hessian  can  describe  the 
dynamically  important  regions  of  the  ground-electronic  state  potential  energy  surface, 
including  the  comer-cutting-tunneling  region  of  the  reaction  swath,  with  practical 
accuracy*^. 

3.4  Error  Estimates 

In  conjunction  with  the  above  referenced  results,  we  also  added  mathematically  based 
error  estimates.  We  used  three  smalt  sets  of  barrier  heights  for  heavy-atom  transfer, 
nucleophilic  substitution,  and  unimolecular  and  association  reactions  as  benchmarks  for 
comparing  and  developing  theoretical  methods.  The  data  sets  were  chosen  to  be 
statistically  representative  subsets  of  the  NHTBH38/04  database.  Each  data  set  consists 
of  6  barrier  heights;  we  called  these  small  benchmark  suites  HATBH6,  NSBH6,  and 
UABH6.  Benchmark  values  were  tabulated  for  204  combinations  of  theory  level  and 
basis  set.  The  theory  levels  studied  includes  single-level  wave  function  theory  like 
Hartree-Fork,  Miller-Plesset  perturbation  theory,  quadratic  configuration  interaction,  and 
coupled  cluster  theory;  they  also  included  multicoefficient  correlation  methods,  local  and 
hybrid  density  functional  theory,  and  semiempirical  molecular  orbital  methods.  The  three 
new  representative  data  sets  were  combined  with  a  previous  representative  data  set  for 
hydrogen-transfer  reactions  to  form  a  new  compact  but  diverse  and  representative  data  set 
called  DBH24. 

Comparison  of  a  large  number  of  methods  for  their  performance  on  DBH24  leads  us  to 
recommend  the  following  methods  for  barrier  height  calculations,  in  order  of  decreasing 
cost:  G3SX,  BMC-CCSD,  PWB6K,  BBIK,  M06-L,  MPWIK,  HF/MIDI!,  and  PM3.  The 
three  small  but  representative  data  sets,  HATBH6,  NSBH6,  and  UABH6,  were  identified 
for  the  barrier  heights  of  heavy-atom  transfer,  nucleophilic  substitution,  and  unimolecular 
and  association  reactions,  respectively.  They  were  representative  of  the  full  data  set 
within  12%  (HAT),  9%  (NS),  and  11%  (UA),  respectively.  We  combined  these  data  sets 
with  a  previous  small  representative  data  set  for  hydrogen-transfer  reactions  to  create  a 
diverse  representative  data  set  of  zero-point-exclusive  barrier  heights  called  DBH24. 
Assessment  of  methods  with  16  DBH24  showed  that  DFT  and  multilevel  methods  have 
much  better  performance-to-cost  ratios  than  single-level  WFT  methods.  The  best  N6 
method  was  BMC-CCSD,  and  its  cost  was  an  order  of  magnitude  smaller  than  the  best 
N7  methods,  although  it  is  almost  as  accurate.  The  two  best  N4  methods,  PWB6K  and 
BBIK,  outperformed  the  best  N5  method  MC3BB.  The  best  local  DFT  method  was  M06- 
L. 


"Optimizing  the  Performance  of  the  Multiconfiguration  Molecular  Mechanics  Method,"  O.  Tishchenko 
and  D.  G.  Truhlar,  Journal  of  Physical  Chemistry  A,  in  press.  (Proofs  returned  to  ACS  Publications  on  Oct. 
31,2006). 
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The  MG3S  basis  set  gave  the  best  performance  for  DFT  methods  among  the  tested  basis 
sets'^. 

4  Modified  POLYRATE  Program:  Test  software 

This  seetion  details  the  final  version  of  the  modified  POLYRATE  software  that  was 
coded  to  implement  the  improvements  described  above.  A  GUI  was  designed  based  upon 
the  user’s  response  to  input  to  the  questions  on  the  GENERAL  INFORMATION  screen 
unless  stated  otherwise.  Additionally,  this  screen  will  also  prompt  the  user  to  enter  the 
number  of  the  reactants  and  products. 

4.1.1  General  Section 

TITLE  Keyword 

This  keyword  is  used  to  give  the  title  for  the  computation  run.  It  will  be  filled  from  the 
user  free-form  input. 

A  TOMS  Keyword 

This  keyword  will  be  filled  in  automatically  from  the  information  the  user  provides  about 
the  reactants  (see  GEOM  keyword  in  0). 

CHECK  Keyword 

This  keyword  specifies  if  the  input  is  checked  for  consistency  prior  to  beginning  the 
calculation.  Current  default  is  OFF,  but  in  the  new  product  this  keyword  will  always  be 
set  to  ON. 

CLASSVIB  Keyword 

This  keyword  specifies  if  vibrational  partition  functions  are  calculated  classically.  It  will 
be  set  to  NOCLASSVIB  (current  default),  indicating  that  all  vibrations  are  treated  as 
quantized. 

DL  Keyword 

This  keyword  determines  if  dual-level  interpolated  corrections  are  done.  It  will  always  be 
set  to  none. 


"Representative  Benchmark  Suites  for  Barrier  Heights  of  Diverse  Reaction  Types  and  Assessment  of 
Electronic  Structure  Methods  for  Thermochemical  Kinetics,"  J.  Zheng,  Y.  Zhao,  and  D.  G.  Truhlar,  Journal 
of  Chemical  Theory  and  Computation,  in  press,  (accepted  Nov.  13, 2006) 
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ICOPT  Keyword 


This  keyword  sets  options  for  interpolated  optimized  corrections.  It  will  not  be  used  in 
the  new  product. 

GEOMUNIT  Keyword 

This  keyword  specifies  the  unit  for  the  geometries.  It  will  be  filled  from  the  user  input  in 
the  pull-down  menu. 

IVTSTO  and  IVTSTl  Keywords 

TBD. 

MDMOVIE  Keyword 

This  keyword  determines  if  the  user  wants  to  use  the  calculation  results  in  an  animation. 
The  default  is  OFF,  but  the  option  to  turn  it  on  will  be  provided. 

RESTART  Keyword 

This  keyword  directs  the  program  to  use  already  existing  files  and/or  save  some  files  for 
future  use.  All  such  options  will  be  available  in  the  new  product. 

SUPERMOL  Keyword 

This  keyword  indicates  if  the  supermolecule  mode  is  used.  It  is  relevant  only  if  the 
potential  energy  is  calculated  from  scratch  (see  below)  and  will  always  be  set  to  default, 
which  is  ON. 

WRITE62,  WRITEFU30,  and  WR1TEFU31  Keywords 

These  keywords  determine  if  certain  files  are  written  out  for  future  use.  All  these  options 
will  be  available  in  the  new  product. 

4.1.2  Energetics  Section 

This  section  will  be  built  from  the  user  input  to  the  questions  in  the  GENERAL 
INFORMATION  screen,  just  like  the  GENERAL  section.  It  is  used  to  specify  the  method 
for  computing  the  potential  energy  surfaces  (PES).  The  choices  are  to  calculate  from 
scratch  using  the  so-called  hooks  or  to  use  already  existing  file  saved  in  previous 
calculations.  These  options  will  be  available  in  the  new  product. 

If  the  user  decides  to  calculate  potential  energy  from  scratch,  the  calculation  will  be 
accomplished  by  a  seamless  interface  with  one  of  the  electronic  structure  packages 
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developed  by  the  University  of  Minnesota,  e.g.,  GAMESSPLUSRATE,  which,  in  turn, 
calls  GAMESSPLUS,  another  package  developed  by  the  University  of  Minnesota. 

It  is  of  importance  to  note  that  no  specialized,  expert-level  knowledge  of  PES  is  needed 
to  provide  input  for  this  section.  Nevertheless,  some  knowledge  of  PES  is  needed  for 
understanding  of  transition  state  theories.  Therefore,  treatment  comparable  to  Chapter  7 
of  the  reference  1  above  will  be  given  in  the  accompanying  tutorial  handbook. 

4. 1 .3  Ancillary  Section 

This  section  contains  information  about  the  method  for  computing  the  Hessian.  It  is 
relevant  only  if  the  PES  calculations  are  performed  from  scratch  using  the  hooks.  If  it  is, 
indeed,  the  case,  this  section  will  be  built  automatically  using  the  current  defaults. 

4. 1 .4  Optimization  Section 

Like  the  Ancillary  Section,  this  section  is  relevant  only  if  the  PES  calculations  are 
performed  from  scratch.  In  such  cases,  an  OPTIMIZATION  screen  will  appear  with 
default  values  pre-filled.  The  user  will  have  an  option  to  override  any  of  the  default 
values.  The  parameters  in  this  section  have  more  to  do  with  numerical  mathematics  than 
with  chemical  kinetics.  In  most  cases  user  will  most  likely  let  the  default  values  stand. 

4. 1 .5  REACT,  PROD,  WELL,  and  START  Sections 

These  sections  describe  the  geometry  and  (optionally)  other  quantum-mechanical 
properties  of  the  reactants  and  the  products.  They  have  the  same  set  of  keywords. 
REACT  and  PROD  sections  must  be  built  for  each  reactant  and  product  respectively,  and 
they  are  named  REACT  1,  REACT2,  PRODl,  and  PROD2.  A  corresponding  number  of 
screens  will  be  provided  for  the  user  to  fill  in  the  information. 

Many  kejwords  in  these  sections  involve  more  advanced  concepts  from  quantum 
chemistry  than  the  user  is  assumed  to  be  familiar  with.  Therefore,  they  will  be  given  a 
detailed  treatment  in  the  book  which  will  accompany  the  software. 

STATUS  Keyword 

This  keyword  must  be  the  first  line  of  the  section.  It  will  always  be  set  to  zero  (default), 
indicating  that  geometry,  energy,  frequencies,  and  eigenvectors  (in  START  section)  must 
be  calculated. 

GEOM  Keyword 

This  keyword  is  used  to  describe  the  geometry  of  the  reactants,  products,  and,  if 
applicable,  saddle  point.  It  is  required  in  all  sections.  Since  STATUS  keyword  is  always 
set  to  zero  (see  above)  the  geometries  will  represent  a  user’s  guess. 
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For  each  reactant,  product  and  saddle  point,  the  user  will  draw  the  structural  formula  with 
the  mouse.  The  program  will  then  guess  the  geometry  assuming  that  the  molecule  is 
planar  and  using  the  covalent  radii  of  the  atoms  involved.  This  information  will  also  be 
used  to  build  the  ATOMS  keyword  in  the  GENERAL  section. 

SPECIES  Keyword 

This  keyword  is  used  to  specify  the  linearity  and  the  phase  of  each  species  as  one  of  the 
following: 

•  Atomic 

•  Linear  reactant  or  product 

•  Nonlinear  reactant  or  product  (default) 

•  Solid-state  reactant  or  product 

•  Linear  (generalized)  transition  state 

•  Nonlinear  (generalized)  transition  state 

•  Solid-state  (generalized)  transition  state 

In  the  new  product  it  will  be  determined  automatically  whether  the  species  is  reactant, 
product,  or  transition  state,  depending  on  which  screen  the  user  is  responding  to.  Other 
choices  will  be  presented  in  a  form  of  the  pull-down  menu,  depending  on  the  entry  for  the 
GEOM  keyword. 

If  the  species  is  single-atom,  the  valid  choices  are  atomic  and  solid-state.  If  the  species  is 
diatom,  then  the  valid  choices  are  linear  and  solid-state.  For  all  types  of  species  valid 
choices  are  linear,  nonlinear  and  solid-state. 

CONSTANT  Keyword 

This  keyword  is  used  to  indicate  coordinates  held  constant  in  optimization.  It  will  not  be 
used  (current  default)  in  the  new  product. 

ELEC  Keyword 

This  keyword  is  used  to  indicate  the  degeneracy  of  the  electronic  states.  Current  default 
(all  states  singly  degenerate)  will  be  used  in  the  new  product. 

ENERGY  Keyword 

This  keyword  is  used  only  if  STATUS  is  set  to  4  or  6.  Since  in  the  new  product  the  status 
will  always  be  set  to  zero,  this  keyword  will  not  be  used. 
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FREQ  Keyword 


This  keyword  indicates  that  normal  mode  analysis  must  be  performed,  which  will  always 
be  done  in  the  new  product. 

FREQUNIT  Keyword 

This  keyword  specifies  user’s  choice  to  express  frequency  in  cm''  or  in  atomic  units.  This 
choice  will  be  given  to  the  user  via  a  pull-down  menu. 

HESSIAN  Keyword 

This  keyword  is  used  only  if  STATUS=4.  Since  in  the  new  product  STATUS  will  always 
be  set  to  zero,  this  ke3word  will  not  be  used. 

INITGEO  Keyword 

This  keyword  will  always  be  set  to  default  value  of  GEOM  indicating  that  the  initial 
geometries  will  be  taken  from  the  GEOM  keyword  (See  above). 

VIB  Keyword 

This  keyword  is  used  only  if  STATUS=6.  Since  in  the  new  product  STATUS  will  always 
be  set  to  zero,  this  keyword  will  not  be  used. 

LINAXIS  Keyword 

This  keyword  is  used  to  show  the  orientation  of  a  linear  molecule.  The  default  value,  z- 
axis,  will  be  used  in  the  new  product. 

PROJECT  Keyword 

This  keyword  is  used  to  indicate  that  projection  operators  are  applied,  which  usually 
leads  to  better  results.  Current  default,  which  is  to  apply  projection  operators,  will  be 
used  in  the  new  product. 

DIA  TOM  Keyword 

This  keyword  is  used  to  treat  diatomic  molecules  with  a  Morse  function.  Current  default 
is  to  treat  diatomic  molecules  in  the  same  way  as  all  the  others.  Since  there  is  a  separate 
set  of  keywords  for  anharmonicity,  it  appears  to  be  almost  redundant.  This  keyword  will 
not  be  used  in  the  new  product. 
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4. 1 .6  Harmonic/ Anharmonicity  Sections  (Optional) 

Several  keywords  are  presently  available  in  POLYRATE  to  indicate  a  choice  to  treat  all 
modes  as  harmonic  oscillators  (HARMONIC  keyword)  or  to  use  one  of  the  following 
anharmonicity  options: 

•  All  modes  Morse  (MORSE  keyword) 

•  All  modes  Morse/quadratic-quartic  (MORSEQQ  keyword) 

•  All  vibrational  modes  are  treated  by  the  WKB  method  with  quadratic-quartic  fit 
(QQWKB  keyword) 

•  All  vibrational  modes  are  treated  by  the  semi-classical  WKB  method  with 
quadratic-quartic  fit  to  potential  (QQSEMI  keyword) 

These  options  will  be  available  to  the  user  via  a  pull-down  menu.  However,  the  user’s 
choice  will  be  treated  as  a  default,  i.e.  there  will  be  an  option  to  treat  some  of  the  modes 
differently. 

For  instance,  ground  state  can  be  treated  by  the  WKB  method  with  true  potential  and 
some  modes  can  be  treated  using  the  hindered-rotor  approximation.  A  user  can  chose  one 
or  both  these  two  options  with  buttons,  in  which  case  additional  pop-up  screens  will 
appear  prompting  the  user  to  fill  in  the  required  information. 

The  following  paragraphs  describe  specific  keywords  used  with  the  anharmonicity 
options. 

ANTLR  Keyword 

This  keyword  specifies  the  criterion  on  the  Morse  anharmonicity  parameter  Xg  and  is  used 
only  with  Morse/quadratic-quartic  calculation.  If  the  user  chooses  this  option,  an 
appropriate  box  will  be  activated  with  the  pre-filled  default  value  of  1x10”*.  The  user 
may  choose  to  override  this  default. 

DEMIN  Keyword 

This  keyword  is  needed  if  any  of  the  Morse  anharmonicity  options  is  chosen.  It  is  used  to 
indicate  the  lowest  dissociation  energy.  Same  as  for  the  ANTLR  keyword,  in  such  cases 
an  appropriate  box  will  be  activated.  POLYRATE  uses  the  default  value  of  0.159  Eh,  but 
it  is  not  recommended.  Therefore,  a  detailed  explanation  of  this  concept  will  be  provided 
in  the  book  that  will  accompany  the  new  product. 

MORMODEL  Keyword 

This  keyword  is  used  to  give  specifics  if  the  user  chooses  Morse  or  Morse/quadratic- 
quartic  model.  There  are  essentially  five  choices  that  can  be  made.  They  will  be  provided 
to  the  user  via  a  pull-down  menu  if  the  choice  is  made  to  use  Morse  or  Morse/quadratic- 
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quartic  model.  Detailed  description  of  these  choices  will  be  provided  in  the 
accompanying  tutorial  handbook. 

4.2  POLYRATE  Tutorial  Handbook 

Dr.  Truhlar  of  U  of  Minn  detailed  the  above  modifications  and  GUI  into  a  POLYRATE 
tutorial  handbook,  explaining  the  algorithms  in  the  current  version  of  POLYRATE.  This 
tutorial  will  be  published  in  Reviews  of  Computational  Chemistry^'^.  A  copy  of  this 
paper,  along  with  all  papers  being  generated  using  and  citing  the  results  of  these  efforts  is 
being  sent  to  the  Government  sponsor  on  CD-ROM.  To  make  POLYRATE  portable,  it 
has  been  modified  to  run  under  LINUX,  specifically,  under  g77  compiler,  the  standard 
compiler  for  LINUX  platforms. 

5  CONCLUSIONS 

This  Phase  I  project  has  developed  modifications  to  POLYRATE  capable  of  being  used 
by  a  chemist  or  a  chemical  engineer  who  is  not  an  expert  in  chemical  kinetics  or  quantum 
chemistry.  The  new  more  complex  software  system  design,  if  implemented,  could  have 
an  even  more  broad-based  user  community.  Through  the  use  of  interactive  GUIs  which 
access  the  newer  higher  order  functionality  of  the  software,  the  system  could  be 
accessible  to  any  user  familiar  with  the  basic  concepts  and  terminology  of  quantum 
chemistry,  statistical  mechanics  and  chemical  kinetics. 

For  the  Phase  2  of  this  effort,  a  seamless,  integrated  computer  code  can  be  developed  to 
permit  the  accurate  calculation  of  rate  constants  for  chemical  reactions  based  on 
electronic  structure  calculations.  The  code  would  be  computationally  efficient  and  would 
contain  a  user  interface  that  permits  easy  use  of  the  program  by  personnel  who  are  not 
experts  in  computational  chemistry.  By  creating  an  easy  to  use,  rate  predictive  code  that 
does  not  rely  on  information  from  existing  databases,  accurate  estimates  for  unknown 
rates  in  kinetic  models  can  be  obtained.  New  materials  can  also  be  able  to  be  more 
quickly  designed  and  synthesized  in  bulk.  The  approach  can  be  useful  for  commercial 
and  research  applications  in  materials,  chemical  vapor  deposition  (CVD),  and  combustion 
chemistry  applications  including  the  automotive,  aerospace,  and  electronics  industries. 


"Variational  Transition  State  Theory  with  Multidimensional  Tunneling,"  A.  Femandez-Ramos,  B.  A. 
Ellingson,  B.  C.  Garrett,  and  D.  G.  Truhlar,  in  Reviews  in  Computational  Chemistry,  Vol.  23,  edited  by  K. 
B.  Lipkowitz  and  T.  R.  Cundari  (Wiley-VCH,  Hoboken,  NJ,  2007),  pp.  125-232.) 


Contract  No.  FA9550-05-C-0065 


18 


